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A trajectory approach is taken to the hydrodynamical 
treatment of collective excitations of a Bose-Einstein conden- 
sate in a harmonic trap. The excitations induced by linear 
deformations of the trap are shown to constitute a broad class 
of solutions that can be fully described by a simple nonlinear 
matrix equation. An exact closed-form expression is obtained 
for the solution describing the mode {n — 0, m — +2} in a 
cylindrically symmetric trap, and the calculated amplitude- 
dependent frequency shift shows good agreement with the ex- 
perimental results of the JILA group. 



The recent realisation of Bose-Einstein condensation 
in dilute trapped atomic vapours Jj]-||] has motivated an 
extensive theoretical study of these systems (for a review 
see reference pj ) . Many of the nonlinear features of the 
condensates are manifested in their collective excitations 
!&-§[. In this letter we consider the collective excita- 
tions of a condensed atomic gas in a harmonic trap in 
the regime in which the number of atoms is sufficiently 
large that the hydrodynamical treatment of Stringari is 
valid H . We apply this treatment by means of a tra- 
jectory approach, which constitutes a generalisation of 
the scaling transformation of Castin and Dum ]l0| . To 
first order in the excitation amplitude our results for the 
oscillatory modes in a stationary trap are identical to 
those of the perturbative treatment of Ohberg et al. [jll| . 
However for a broad class of excitations, our analysis pro- 
duces a simple nonlinear matrix equation, which is valid 
to all orders in the excitation amplitude. The excitations 
to which this matrix equation applies are those induced 
by time-varying linear deformations of the trap (that is, 
perturbations of the trapping potential that maintain its 
harmonicity) . This class includes all the excitations that 
have to date been studied experimentally. 

We suppose that the condensed gas is confined by a 
harmonic potential, which, in its most general form, is 
given by 

U(r, t) = lj2 M*) N ~ n(t)] [rj ~ fj(t)] , (1) 

where f(t) is the position of the centre of the trap, 
and JC(t) is a symmetric matrix with components /Cy(t), 
which represent the spring 'constants' of the trap. 



In the mean field approximation, the evolution of 
the condensate wavefunction is determined by the time- 
dependent Gross-Pitaevski equation, which, for zero tem- 
perature, is 



ihd t $ (r,i) = 



-—A + U(r,t)+Ng\<S>(r,t)\ 2 



$(r,t) 
(2) 



where N is the number of atoms in the condensate, and 
g is a measure of the strength of the atomic interactions. 
g is related to the s-wave scattering length a through 
g = 4:nh 2 a/m, and is assumed to be positive. 

Defining the density of the condensate as p (r, t) — 
N |$ (r,£)| 2 , and its velocity field as v(r,t) = 
(h/2im) Vlog[$(r,t)/$*(r,t)], it is easily shown that the 
Gross-Pitaevski equation (^) is equivalent to the follow- 
ing equations for p and v 



dtp + V (vp) = 



1 



md t v + V ( U eS + ^mv 2 ) = 0, 



where 



U c s = U + gp 



h 2 



Im^fp 



Av/p 



(3) 
(4) 

(5) 



is an effective potential. 

Following Stringari |J we consider the limit in which 
the number of atoms N is sufficiently large that the den- 
sity p becomes smooth, and it is valid to neglect the 
kinetic energy pressure term h 2 (Ay/p) jlra^fp. The ra- 
tio between the typical force caused by this force and 
the force — V(gp) caused by inter-particle interactions 
can be shown to be of the order of {h 2 jmR^)lp, ~ 
N~ 2 / 3 (T C / p,) 2 , and for the most of the BEC experiments 
this ratio is as small as 1CP 1 — 1CP 3 . (Here and below 
Rq ~ \J p/muj 2 is the Thomas-Fermi radius of the con- 
densate, T c ~ huN 1 / 3 is the Bose-Einstein condensation 
temperature, p ~ gp is the chemical potential, and to is 
a typical trapping frequency.) In this regime the con- 
densate can be modeled as a classical gas in which each 
particle is subject to a force 



F(r,i) = -V r [U{r,t)+gp{r,t)] 



(6) 
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and for which the velocity field is irrotational. We include 
the subscript r on the gradient operator V to indicate 
explicitly that all derivatives are taken with respect to 
the components of r. This notation is used since we shall 
later introduce a change of coordinates. 

The requirement that V r xv(t) = imposes a con- 
straint only on the initial condition, the reason being 
that since the only force involved (Q) is the gradient of 
a potential, the total derivative of V r x v(<) vanishes. 
This means that, provided the velocity field is initially 
irrotational everywhere, it will remain so with time. 

Using expression (||) for the classical force, we obtain 
the following evolution equation for r(t) 



mf (t) = -Kit) [r(t) - r(t)] - gV r p (r, t) . 



(7) 



In order to solve this equation, it is necessary to deter- 
mine the density of the condensate at position r and time 
t, which in turn requires a knowledge of the trajectories 
of all the other particles in the condensate. To this end, 
we associate with each trajectory r{t) a unique reference 
point i"o , corresponding to the position of a particle in the 
ground state $o of the condensate in the hydrodynami- 
cal limit for a stationary harmonic trap Uq(t), given by 
equation (Q) with r(i) = and K.(t) = Kq. $o is obtained 
from the solution p of equation (g) with F(r,i) = 0, 
through 



/9o(r ) = iV|* (r )r 



and is given by 



$o (r ) 



Ifi - U (r ) 



Ng 



(8) 



(9) 



where \i is the chemical potential. Note that this is just 
the Thomas-Fermi solution. 

The density of the condensate at time t can now be 
written as 



p(r,t) = dct 



dr 
dr 



Po (r ) ■ 



(10) 



where dr /dr is a matrix whose ij th element is dr, / dr j . 



Given that V r = 



(dr/dr y 



V r „ , we can rewrite 



equation (0) as a partial differential equation for r(ro,i) 
md 2 t r = -K(t) [r - f(t)} 
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dr 
dr 



del ( ) /.nU'ni 



(11) 



The functions r satisfying this equation can be used 
to construct hydrodynamical solutions of the Gross- 
Pitaevski equation @. Given that the velocity field is 



irrotational, we can express it as the gradient of some 
function (r, t) 



v(t) = V r 0(r,f) 
in terms of which we obtain 



$(r,i) 



e -i/3(t) e ime(r,t)/h 



$o (r ) 



^/det (dr/dro) ' 



(12) 



(13) 



where ro is uniquely defined in terms of r and t, and (3(t) 
is some function of time. 

By applying this analysis to oscillatory modes in the 
stationary harmonic potential U(r,t) — Uo(r) (that is, 
K(t) = ICq, r(t) = 0) in the limit of small excitation 
amplitudes, we recover the results of Ohberg et al. | jll| , 
as we now show. For low amplitude oscillations 

r (r , t) « r + e [St (r ) e^ 1 + 5r* (r ) e^] , (14) 

where e is a perturbation parameter. Substituting this 
expression into equation ( pd| ) , and keeping only first order 
terms in e we obtain the relation 

muj 2 5r = V {8r ■ (/C r) — [// — U (r)] V ■ 5r} , (15) 

in which we have dropped the distinction between r and 
ro, since this affects only higher orders. 5r is clearly 
expressible as a gradient 

5r = VW, (16) 

where the function W satisfies the equation 

muj 2 W = VW ■ (K t) U (r)} AW. (17) 

Expanding the solution ( |l3|) to first order, we obtain 



$(r,t) 
where 

and 



[$o (r) + e (we" 



U + f- 



/+-/- 



)], (18) 



(19) 



1 - 



C^o(r) 



±1/2 



w. 



(20) 



with C-/C+ = fuo/2u. The expression (|T^), with the 
differential equation (|7|) for W, is precisely the result 
obtained by Ohberg et al. [^T). We have thus shown 
that their solution describes the elementary phonon-likc 
excitations of the nonlinear hydrodynamical model. 

The oscillatory modes for which the function W(r) is 
of third or higher order in the components of r can- 
not be excited by linear deformations of the harmonic 
trap. That is, they cannot be produced by applying a 
time-dependent harmonic potential of the form (Q) to 
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the ground state condensate $o ( r ) ■ However, the modes 
for which W(r) is linear or quadratic can be produced 
in this way, and for these modes it is possible to go be- 
yond the perturbative limit, to obtain a simple nonlinear 
matrix equation that fully describes the dynamics of the 
condensate in the hydrodynamical regime. 

If the condensate is initially in the ground state $o (r), 
then the effect of a time-dependent harmonic potential 
(|l|) is simply to distort the condensate in a linear fashion 



r(r 0> t) =R(t) + M(t)r . 



(21) 



Here R(t) is the centre of mass of the condensate, and 
M it) describes the time-dependent contraction or shear- 
ing of the condensate along various directions. The ma- 
trix dr/dro is now independent of position 



dv 

— — = M(t) = constant (ro), 

OTo 



(22) 



and equation ( |ll|) can be separated into a part that is 
independent of position, and a part that depends linearly 
on position, giving the following two relations 

mR = —K(t) [R — r(i)] (23) 
mM = —K{t)M + dot (My 1 (M~ 1 ) T /C . (24) 

From equation (]23| ) it is clear that for a stationary trap 
(lC(t) = /Co) the centre of mass motion is separable along 
the three axial directions, and the frequencies of oscilla- 
tion along these directions are independent of amplitude, 
in agreement with the Kohn Theorem (see pl| and ref- 
erences therein). 

The velocity can be written in terms of r as 



v = R + 7W7W- 1 (r-R), 



(25) 



and hence the requirement that the velocity field have 
zero curl is equivalent to the constraint that the matrix 
MM^ 1 be symmetric 



MMT 



(26) 



Using equation ( p4|) it is easy to verify that provided 
condition ( p6|) is satisfied at the initial time it will be 
satisfied at all later times. 

Given this constraint it is clear that for a stationary 
potential in the limit of small excitations, there must 
be exactly six linearly independent modes for which the 
function W(r) is quadratic in the components of r. Fol- 
ey lindrically symmetric traps (lOx = ^Oy = they are 
(in the notation of reference {n = 2, m = 0, (±)}, 

{n = 0, m = ±2} and {n = 1, m = ±1} |l3|. In this small 
amplitude limit we can write 



M(t) = l + e[SMe 



8 Me 1 ' 



(27) 



where SM is a constant matrix, independent of both posi- 
tion and time. It completely characterises the modes and 
their linear combinations in the small-amplitude limit. 
For the scaling modes described by Castin and Dum |Tc| ] 
(the modes {n = 2, m = 0, (±)}, and a symmetric super- 
position of the modes {n — 0, m = ±2}) 5M is diagonal. 
For the mode {n = Q,m = +2}, which we shall consider 
next, it is 



SM 




(28) 



An exact solution for this mode can be found by setting 
IC(t) = /Co in equation (p4|), and is given by 



M(t) =exp [\Q(t)]K x (t), 



where 



C(*) = 



cos (ujt) 
sin (ut) 




sin (ojt) 
- cos (ujt) 




(29) 



(30) 



and lo is the mode frequency. From expression ( J3 0| ) for 
Q(t) it can be seen that the transformation exp [AQ(t)] 
describes the distortion of the condensate into an ellipse, 
and the rotation of that ellipse about the z axis. The 
parameter A is a measure of the amplitude of excitation, 
and is precisely the quantity we labelled e in the small 
amplitude limit. 

The matrix lZ\(t) describes a slow rotation about the 
z axis, together with a centrifugal stretching in the x-y 
plane, and a corresponding contraction along z 



/3cos (uj s t) 
(3 sin (ui s t) 




— /3sin (u s t) 
(3 cos (ui s t) 







-2/3 



(31) 



This rotation is required in order for M (t) to satisfy con- 
dition ( p6| ) . 1Z\ (t) depends on A through both the param- 
eter (3 and the slow frequency uj s . 

The aspect ratio a of the condensate is defined as 



Li — L s 
U + L s 



(32) 



where Li/L s is the ratio of the lengths of the long and 
short axes of the condensate in the x-y plane. It is 
related to the parameter A by a — tanhA. The con- 
stant P depends in turn on a through f3 = [(1 + 



a 4 )/(l 



The slow frequency co s is given by 



lu s = [\/2a 2 /Vl + a 4 ]^_L, and the mode frequency lo 
scales with the aspect ratio as 



V2- 



zLO±. 



(33) 



lo thus depends on the amplitude of excitation. For very 
low amplitudes, A — > 0, a — > 0, and the mode frequency 
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tends to uj — > \/2ui±. At very high amplitudes, A — > oo, 
a — > 1 and the mode frequency approaches the limit w — > 
2ui_ , which coincides with the frequency of oscillation of 
a non-interacting gas. 

This mode has been observed experimentally by Jin 
et al. HD, wno measured its frequency as a function of 
response amplitude (aspect ratio) after opening the trap 
and allowing the condensate to expand. Due to the pres- 
ence of atomic interactions, the aspect ratio will change 
during the expansion. In order to compare our results 
with the experimental data it is therefore necessary to 
model the evolution of the condensate after the trap has 
been opened. We assume that the opening of the trap 
occurs instantaneously. Taking the solution ( ]29| ) for the 
trapped condensate as an initial condition, we insert it 
into (|24|), setting K{t) = 0, and integrate numerically. 
The aspect ratio of the condensate is calculated from the 
elements of M. after the appropriate expansion time. 




20 40 60 

response amplitude [%] 



FIG. 1. The normalised mode frequency w/w(0) is plotted 
against response amplitude (aspect ratio a) after expansion. 
Our theoretical predictions (curve) are compared with the 
experimental data of Jin et al. 

Shown in Figure 1 is the experimental data obtained 
by Jin et al. Q for the mode frequency as a function of 
aspect ratio. These points were obtained with a radial 
trap frequency uj±/2it of 132 Hz and an expansion time 
of 7 ms ||, and the frequencies were normalised to the 
small amplitude value w(0) |l4|]. The ratio of the axial to 
radial trap frequencies was uj z /lli± = \/8. The number of 
condensed atoms was approximately 4500, and temper- 
ature was T = (0.50 ± 0.08)T C , where T c is the critical 
temperature. The theoretical curve, calculated for the 
same parameters, shows good agreement with the exper- 
imental data. 

Note that in our treatment we totally neglected the the 
presence of the non-condensed particles at a finite tem- 
perature T. Within the mean-held framework the force 
caused by the non-condensed particles can be shown to 
be (T/T c ) 3/2 (fi/T c ) 3/2 times smaller than the one related 
to the condensate, which makes this assumption valid up 



to a fraction of the transition temperature. 

In conclusion, we have considered the collective excita- 
tions of a Bose-Einstein condensate in the hydrodynam- 
ical regime. We have shown that the excitations induced 
by arbitrary time-dependent linear deformations of the 
trapping potential form a broad class of solutions of the 
hydrodynamical equations, which remains closed under 
evolution. The dynamics of the condensate for this class 
of excitations can be described by a nonlinear equation 
involving a 3 x 3 matrix. We have obtained an exact 
solution of this equation for the {n = 0, m = +2} mode, 
which shows good agreement with experimental results. 

We would also like to note that described in the liter- 
ature 3-fold class of nonlinear solutions Jl(| covers two 
{n — 2, m = 0} modes and any linear combination of 
{n = 0, m = ±2} modes with (generally complex) equal 
by the absolute value coefficients: the observed in the 
experiment || {n = 0, m = +2} mode does not belong 
to this class. The 6-fold closed class of nonlinear solu- 
tions obtained in our work extends the treatment to two 
other modes {{n — 0, m = ±2} respectively) as well as 
to an arbitrary superposition of {n = 0, m = ±2} modes, 
including {n = 0, m = +2}. 

In the other words, the existing treatment can de- 
scribe any excitation obtained by a modulation of the 
matrix elements of the spring-constant tensor ICij (t) (see 
(|l|)) provided it always remains diagonal. This model can 
not be applied to describe the rotation-like modulations 
such as the "m = 2" JILA experiment ||. Instead, our 
paper presents an extention for this model which can be 
applied to any kind of the time dependence of the spring- 
constant tensor, including the one used in the above ex- 
periment. 
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